%%%%%%%%%% Function FermiFunction %%%%%%%%%%
    function F = FermiFunction(E,Ef, kT)
    he = E(2)-E(1);
    Ne = length(E);
    Nx = length(Ef);
    
    % kT ~= 0, a smoothened distribution
    if kT ~=0
        % Ef is a vector (e.g. Ef + V(x) )
        if length(Ef)>1
            E_extend = E'*ones(1,Nx);
            Ef_extend = ones(Ne,1)*Ef;
            F = 1./(1 + exp((E_extend-Ef_extend)/kT));
        % Ef is a scalar
        else
            F = 1./(1 + exp((E-Ef)/kT));
        end
    % kT =0, a step function
    else 
        % Ef is a scalar
        if length(Ef)>1
            F = zeros(Ne,Nx);
            index = round((Ef-E(1))/he);
            F(1:index,1:end) = 1; 
        % Ef is a vector (e.g. Ef + V(x) )
        else
            index = round((Ef-E(1))/he);
            F = zeros(Ne,1);
            F(1:index) = 1;
        end
            
    end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%